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ABSTRACT 

Aims. We present the results of global 3-D MHD simulations of stratified and turbulent protoplanetary disc models. The aim of this work is to 
develop thin disc models capable of sustaining turbulence for long run times, which can be used for on-going studies of planet formation in 
turbulent discs. 

Methods. The results are obtained using two codes written in spherical coordinates: GLOBAL and NIRVANA. Both are time-explicit and use 
finite differences along with the Constrained Transport algorithm to evolve the equations of MHD. 

Results. In the presence of a weak toroidal magnetic field, a thin protoplanetary disc in hydrostatic equilibrium is destabilised by the mag- 
netorotational instability (MRI). When the resolution is large enough (~ 25 vertical grid cells per scale height), the entire disc settles into a 
turbulent quasi steady-state after about 300 orbits. Angular momentum is transported outward such that the standard a parameter is roughly 
4 - 6 X 10"'. We find that the initial toroidal flux is expelled from the disc midplane and that the disc behaves essentially as a quasi-zero net 
flux disc for the remainder of the simulation. As in previous studies, the disc develops a dual structure composed of an MRI-driven turbulent 
core around its midplane, and a magnetised corona stable to the MRI near its surface. By varying disc parameters and boundary conditions, we 
show that these basic properties of the models are robust. 

Conclusions. The high resolution disc models we present in this paper achieve a quasi-steady state and sustain turbulence for hundreds of 
orbits. As such, they are ideally suited to the study of outstanding problems in planet formation such as disc-planet interactions and dust 
dynamics. 
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1. Introduction 

Observational surveys of star forming regions in the Galaxy 



and dust orbiting i 


/oung stars (e.g 


. .Beckwith & Sargent 199^; 


O'dell et all 


1 19931 


IStauffer et alJ 


1 19941 ISicilia-Asuilar et alJ 


2006VKessler-Silacci et al.'2006V 


It is commonly believed that 



Il969t E issaue 3 ll993l) . The discovery of numerous extrasolar 
planets has increased the need for greater understanding of their 
properties so that accurate models of planet formation can be 
developed. 

These "protoplanetary" discs often show evidence for ac- 
tive accretion with a canonica l mass flow rate onto the ce ntral 
star of ~ 10"** Mo yr"' (e.g. ISicilia-Aguilar et aDl2004 . re- 
quiring a source of anomalous viscosity to transport angular 
momentum outward. It has long been beli eved that this is pro- 
vided by turbulence within the disc (e.g. IShakura & SunvaevI 
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[19731). So far only one mechanism has been shown to work 
reliably: MHD turbulence generated by the magnetorotational 
instab ility (MRI) (Balbus & Hawley 1 991; Hawl ev & B albua 
Il99ll) . Given the cool and dense nature of protoplanetary 
discs, there are questions about the global applicability of the 
MRI in such en vironments as the ionisation fraction is low 
jBlaes & BaibusI Il994ll) . Models suggest that protoplanetary 
discs are likely to have both magnetically active zones, where 
the disc is turbulent, and adja cent magnetical l y 'dead-zones' 
where the flow is laminar (e .g. lGammidll996t iFromang et alJ 
l2002t lllgner & Nelso n"200^. In this paper we focus on ideal 
MHD simulations of protoplanetary discs. We will examine the 
dynamics of 'dead-zones' in future work. 

Non linear numerical simulatio ns performed using the 
local shearing box formalism (e.g. 'Hawl ev & BaibusI 119911: 
Hawl ev etalJl 19961 11 Brandenburg et al. 1996:D have shown that 
the saturated non linear outcome of the MRI is MHD turbu- 
lence having an effective viscous stress parameter a between 
~ 5 X 10"^ and ~ 0.1, depending on the magnetic field configu- 
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ration. Outward angular momentum transport can thus occur at 
the rate required to match ob served accretion signat ures onto T 
Tauri stars (see, for example. iHartmann et al .119981 who quote 
a ~ 0.01 as being suggested by the observations). Much of this 
early simulation work was performed in discs with no vertical 
stratification, and so was useful in determining the nonlinear 
outcome of the MRI, but did not provide insights into the global 
structure of these discs in either the radial or vertical direc- 
tion. The question of their vertical structure was addressed by 
^tone et al. (1996) and Miller & Stone (2000) who performed 
shearing box simulations of vertically stratified discs. A ba- 
sic result to come out of these studies is that the discs evolve 
to a structure consisting of a dense region around the miplane 
where turbulence is driven by the MRI, sandwiched between 
a tenuous and magnetically dominated corona which is highly 
dynamic, but stable against the MRI. 

Gl obal MHD simulations of turbulent discs (e.g.lArrpitage l 

199?, 'Hawlev' '2000'> '2001] LSteinacker & Papaloizod l2002^ 
Papaloizou & Nelson 2003) confirm the basic picture provided 
by the local shearing box simulations. These early global sim- 
ulations considered the radial structure of turbulent disc mod- 
els, but employed non stratified cylindrical disc models. Recent 
work has been performed examining the dynamics of global, 
vertically stratified turbulent disc models, and has focussed 
on thick accretion tori aro und black hole s , using either the 
Paczynski-Witta potential llHawlevI EoOflt iHawlev & KroliS 
1200 U iHawlev etalJ l200l') or simulating accretion flows in 
the full Kerr metric (e.g. De Villiers et al. 2003). The starting 
conditions for these models are usually thick, constant angu- 
lar momentum tori for which the gas is assumed to be non- 
radiating. These quickly evolve into thick accretion discs for 
which H/R > 0.1. To date, there have been no published sim- 
ulations of vertically stratified thin discs undergoing MHD tur- 
bulence (i.e. discs more akin to protoplanetary discs). In part 
this is because of the increased computational burden associ- 
ated with simulating thin discs due to the higher resolution re- 
quirements. 

The aim of this paper is to present and analyse a suite of 
MHD simulations of global, stratified, turbulent, protoplane- 
tary disc models with H/R < 0.1. We focus on the global struc- 
ture of the discs, and analyse how modifications to the bound- 
ary conditions and disc thickness affect the results. The longer 
term goal is to use these models to examine a range of prob- 
lems in planetary formation in the context of turbulent, verti- 
cally stratified discs. As such this paper is the first in a series. 
Future publications will address problems such as: the evolu- 
tion of dust in turbulent discs; the orbital evolution of plan- 
etesimals and low mass planets; gap formation and gas accre- 
tion by giant protoplanets; the dynamical evolution of dead- 
zones. A number of previous studies have addresse d these is- 
sues using cylindrical models: Nelson & Paoaloizou (2003); 
Papa loizou et al. (2004); Nelson & Papaloizou (2004); Nelson 
( 1200 5 ) studied gap formation and planet formation in turbu- 
lent discs. They showed in particular that type I migration be- 
comes stoch astic because of the turbule nt density fluctuations 
in the disc. 'Fro mang & NelsonI ( l2005h also used cylindrical 
models to study the radial migration of solid bodies due to 
gas drag. They found rapid accumulation of meter size bod- 



ies in anticyclonic vortices apparently resulting from the tur- 
bulence. Another key problem in planet formation, dust set- 
tling in turbulent discs, has been studied recently using shear- 
ing box models (Johansen et al. 2006; Fromang & Papaloizoiy 
i2006HTurner et al.l2006il) . Because of the local approach, these 
analyses had to ignore radial drift. All of these issues are likely 
to be affected by the simultaneous treatment of radial and ver- 
tical stratification that we consider in this present work. 

The plan of the paper is as follows: in section|2 we present 
the basic equations, notations, and diagnostics we use in this 
work. The set-up of our simulations is detailed in section|3]and 
we present their results in section 0] We focus in particular on 
their global structure, and sensitivity to numerical issues such 
as resolution and boundary conditions. Finally, in section|5] we 
summarise our results and highlight future improvements that 
will be added to the models. 

2. Basic equations 

The equations we seek to solve are the standard equations of 
MHD in a frame rotating with the angular velocity Qrot- In 
Gaussian units these may be expressed as: 

dp 



dt 



+ V.(pv) = 0, (1) 



dv 

— + (VV)V -I- 2£lrotXV 

ot 



--VP- VO-H — (VxB)xB, (2) 

p Anp 

dB 

— = Vx(vxB - 77VxB) . (3) 

at 

The symbols have their usual meanings: p is the density, v is 
the velocity, P is the gas pressure, B is the magnetic field, rj 
is the magnetic diffusivity, and is the sum of the gravita- 
tional and centrifugal potentials. In the following discussion 
we will use two systems of coordinates: cylindrical coordinates 
{R,(f),Z) and spherical coordinates {r,0,(p). Using the later, O 
becomes 



GM 



(4) 



where G is the gravitational constant and M the mass of the 
central protostar To complete equations Q-Q, we use a lo- 
cally isothermal equation of state: 



P = ciRYp 



(5) 



where c(R) is the sound speed which is a fixed function of po- 
sition. 

2.1. Averaged quantities and connection to viscous 
disc theory 

When presenting the results of our calculations we make use of 
a number of averaged quantities. For a quantity Q, we define 
the average Q{R, t) through 



Q(R,t) 



J J pQ(R,(p,Z,t)dzd(p 
-^JJpdzd(f> 



(6) 



where Acp - ^„ax - (pmin is the size of the azimuthal domain 
(which is less than 2n in the disc models we present later in 
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this paper). The integrals are taken over the total disc height 
and azimuth. We define the disc surface density through 



2W 



pdzd(p 



(7) 



Note that without loss of generality we can also include a time- 
average within our averaging procedure. In this case the quan- 
tity Q{R, t) becomes: 



2 r^l+ht 
I I 

2Af J,-A/ 



Q(R,t)dt 



where 2Af is the time interval over which the averaging proce- 
dure is performed. 

We now consider the connection between the tur- 
bulen t disc models and sta n dard viscous disc theory 
(e.g . Shakura & Sunvaev 'l 973t iBalbus & Papaloizoul 1 19991 
IPapaloizou & Nelsoa.2003i) . Averaging the continuity equa- 
tion ([ij gives 



dt RdR 



(9) 



Consider the azimuthal component of the momentum equa- 
tion Multiplying by R to give an equation for the conser- 
vation of angular momentum about the Z axis, and averaging 
over the disc height and azimuth gives: 



RdR 



R^I.\ 



(10) 



Here j is the specific angular momentum, Rv^, and we have 
used the relations vr - vr + 6vr and = -i- (Jv^, where 
6vr and dv^j, are the fluctuations in the radial and azimuthal ve- 
locities. Note that we have neglected terms due to the pressure 
gradient in equation dlOt as the transport of angular momen- 
tum within the disc is primarily due to Reynolds and Maxwell 
stresses. These are given by the two terms on the right hand 
side of equation (I10> . and we define the Reynolds stress by 



Tr = 1L6v^6vr 

and the Maxwell stress by 



'B^Br 



(11) 



(12) 



We define the viscous stress parameter a through 

Tr — T M 

a - — 

P 



where P is the averaged gas pressure. In the discussion of our 
simulation results presented later in this paper we refer to am 
and ur, which are the Maxwell and Reynolds stress contribu- 
tions to the total a value defined in equation ( I13t . We also dis- 
cuss the volume averages of a, and or which we denote as 
(a), (qm) and (aR). 



Combining equations (|9j and (I10> gives: 



'5f 



^dR 



1 d 

2^:^ H- Svr-^ = — — (rWr - R^Tm] 



(14) 



Assuming that the averaged specific angular momentum jiR) is 
time independent, we obtain an equation for the averaged radial 
velocity that is the equivalent of a similar expression obtained 
in standard viscous thin disc theory: 



(8) -.R-^ 



(dT 

BR 



dR^ 



i R- K i M 



)• 



(15) 



By obtaining time averaged values of Tr(R), Tm{R), j(R), 
vr{R) and 2(/?) from our numerical simulations, we are able to 
compare the value of vr(R) obtained directly in the simulations 
to that predicted by equation il5\ . Hence we can examine how 
well our 3-D turbulent and stratified protoplanetary disc mod- 
els are described by thin disc theory, subject to an appropriate 
time average. When calculating these averages in the simula- 
tions, which are performed using spherical polar coordinates, 
we perform the vertical integration by replacing liz — > r sin OdO, 
and integrate along the meridian at fixed r and (p. 

3. Numerical Simulations 

We used two codes t o solve the MHD equati ons described in 
section H GLOBAL ("Hawle v & Ston3ll995 ) and NIRVANA 
fziesler & Yorke 1997). GLOBAL, originally written in cylin- 
drical geometry, was modified to operate using spherical coor- 
dinates. Both codes are time explicit, use finite-differences to 
calculate spatial derivatives and the Method of Characteristics 
Constrained Transport (MOCCT) algorithm to evolve the mag- 
netic field. They have been widely tested in the past on a variety 
of different problems relating to turbulent protoplanetary discs, 
making them particularly well suited to undertake the computa- 
tionally demanding problem of simulating stratified protoplan- 
etary discs models. 



3.1. Stratified disc model set-up 

At time f = we specify a spatial distribution for the hy- 
drodynamic variables that is as close as possible to a strati- 
fied thin disc in hydrostatic equilibrium. Except for the vertical 
stratification, the model pro perties are very close to those of 
IPapaloizou & NelsonI ( l2003h . The mass density p and angular 
velocity Q are defined using: 



(Ro\ 



3/2 / ^2 



(13) n = 



(16) 



(17) 



The radial and meridional velocities and vg are given small 
random perturbations (using a uniform distribution with ampli- 
tude 2.5 % of the sound speed). The disc semi-thickness H is 
related to disc parameters through 



H = 



c(R) 



(18) 
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Model 


Co 


Resolution 


Azimutal extent 


Vertical extent 


Inner radial BC 


Outer radial BC 


Vertical BC 


Code 


Sla 


0.07 


(272, 90, 126) 


;r/4 


0.3 


Reflecting 


Reflecting 


Outflow 


GLOBAL 


Sib 


0.07 


(272, 90, 126) 


nIA 


0.3 


Reflecting 


Reflecting 


Outflow 


NIRAVANA 


S2 


0.07 


(455, 150,213) 


71/4 


0.3 


Reflecting 


Reflecting 


Outflow 


GLOBAL 


S3 


0.07 


(455,150,285) 


7T/4 


0.4 


Reflecting 


Reflecting 


Outflow 


GLOBAL 


S4 


0.07 


(456,150,210) 


nIA 


0.3 


Reflecting 


Reflecting 


Periodic 


NIRVANA 


S5 


0.1 


(360,120,210) 


71/4 


0.43 


Outflow 


Reflecting 


Outflow 


NIRVANA 


CI 


0.07 


(455,150,40) 


7T/4 


0.42 


Reflecting 


Reflecting 


Periodic 


GLOBAL 


C2 


0.07 


(260, 152, 44) 


7T/2 


0.28 


Reflecting 


Reflecting 


Periodic 


GLOBAL 



Table 1. Properties of the models described in this paper. The first column gives the model label, the second gives the sound 
speed at r - tq. The third, fourth and fifth columns describe the resolution and the extent of the numerical box. The type of 
boundary conditions we used are given in columns 6 ,7 and 8, while the last column indicates the code used to run that particular 
model. For the detailed set-up of model C2, see section lT^ or Fromang & Nelson (2005) . 



In the absence of magnetic fields, the initial disc model is com- 
pletely determined once the function c{R) is given. In the fol- 
lowing, we used: 

c(/;) = co(y) (19) 

such that H/R is constant in the models. Using these relations, 
the surface density S satifies 

where So = IpqcqRI'^/ ^^JGM. 

The constants appearing in the above equations are given 
values GM — \, Rq - 1, po - 1 ™d Q^, = 0.5. Depending 
on the models, cq is either equal to 0.07 or 0.1 (see table Q. 
The computational domain extends from /?,„ = 1 to Rout = 8 in 
the radial direction. The vertical and azimuthal extent depend 
on the model. But in all of them, at least 8.5 scale heights are 
covered in 6 in order to provide a good description of both the 
disc midplane and corona. When discussing simulation results 
in this paper, time will be quoted in units of the orbital time at 
the inner edge of the computational domain. With our defini- 
tions, a time span of 500 orbits at r = /?,„ (which is the typical 
duration of a given model) corresponds to 177 orbits at r - 2, 
63 orbits at r = 4 and 27 orbits at r = 7. 

Before adding the magnetic field, the above model 
was run as a purely hydrodynamic disc. Using reflecting 
boundary conditions in r, periodic boundary conditions in 
(f> and outflow boundary conditions in 6, we found it to 
stay very close to the initial model description presented 
above. In particular, we found no sign of the transient 
growth of hy drodynamic perturbations recently proposed in 
the literature (loannou & Kakouris 200 1 : " aF shordi et al]l2005t 
iMukhopadhvav et al.-200S) . The kinetic energy then decays on 
a time scale of ~ 10 orbits (some wavelike motions are then 
excited by the imperfect nature of the boundary conditions on 
time scales of hundreds of orbits, but the associated kinetic en- 
ergy always remains smaller than in the MHD case by at least 
an order of magnitude). 

In the MHD runs, a toroidal magnetic field was added to the 
disc at f = 0. It is defined to be nonzero in the region of the disc 



satisfying 2.5 < r < 6 and \0 - 7t/2\ < C?- We took 6^"/ = 0.1 
in all of the models, except for model S5 for which (^"f = 0.2. 
The strength of the field is such that the ratio p of the thermal 
pressure to the magnetic pressure is everywhere equal to /5o - 
25. Previous global models of such configurations have proved 
to be unstable to the MRI and to generate MHD turbulence that 
spreads into the entire computational domain. We emphasize 
here that the flux of this magnetic field in the (p direction is 
nonzero at the beginning of the calculation. 

3.1.1. Boundary conditions 

We now come onto the boundary conditions. In principle it is 
possible to use reflecting or outflow boundary conditions at the 
inner and outer disc radii. Reflecting boundary conditions are 
unsuitable as waves excited by the turbulence then tend to rat- 
tle around the computational domain in an unphysical manner. 
Outflow conditions were tried in test calculations but proved to 
be unsuitable because of excessive mass loss out of the compu- 
tational domain. Instead we added a non turbulent buff'er zone 
close to each radial boundary (for all runs except S5, see be- 
low), the inner one extending from r = Ri„ to r - 2, the 
outer one extending from r = 7 to r = Rout- In both buffer 
zones, we included a linear viscosity q''" with coeflicient Ci 
(see the definitions of Stone & Norman 1992), to damp fluid 
motion, and resistivity, 77, to create a region that is stable to the 
MRI and therefore non turbulent near the boundary. Both dis- 
sipation coefficients increase linearly from the boundary of the 
buffer zones toward the boundary of the computational domain. 
In all cases except model S3, their maximum value at the inner 
boundary is Ci = Ci,,„ - 1 and rj - r]i„ - 4.9 x 10 while they 
reach Ci = Ci_o«/ - 5 and 77 = rjoui = 10"^ at the outer bound- 
ary. In model S3, we used Cij„ = 1, C^out = 30, rii„ -2x 1Q~'^ 
and rjout = 10"^. Using this set-up, we found that the turbulent 
velocity fluctuations and magnetic stresses damp smoothly in 
the buffer zones. 

For very long integration times, we found that mass tends 
to accumulate at the interface between the inner buffer zone 
and the active disc. This is expected for such a closed bound- 
ary condition. To try to solve this problem, we modified the 
inner boundary condition in model S5 and used a 'viscous out- 
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flow condition' (the outer radial boundary is kept unchanged). 
The viscous outflow condition specifies the radial outflow ve- 
locity in the inner radial boundary using the expression = 
-3/2ac^/Q. there. A value of a = 5 x 10"^ was adopted, in ba- 
sic agreement with the transport coefficient resulting from the 
nonlinear evolution of the MRI (see below). The magnetic field 
boundary condition at the disc inner edge defines the field to be 
normal to the boundary with magnitude defined by the V.B = 
condition (e.g. lHawlevl200ol) . 

The boundary conditions in the 6 direction are also of 
some importance in constructing stratified mode ls. In local 
simulations using the shearing box approximation. ^ Stone et al. l 
( 1996) reported very little effect of the boundary conditions. 
Nevertheless, we investigated their importance by running two 
types of boundary conditions: outflow and periodic. The latter 
is less physical but has the advantage of preserving the total 
flux of the magnetic field and the vanishing value of its diver- 
gence. The former was implemented in two different ways: in 
GLOBAL, a zero gradient boundary condition was applied on 
all th e variables, including the magnetic field ( iMiller & Stond 
I2OO O). NIRVANA, on the other hand, forces the magnetic field 
to become norm al to the bou ndary and still satisfy the condi- 
tion V-B = jHawlevlbOOfll) . The comparison between these 
different alternatives gives an insight to their importance on the 
global structure of the disc. 

In the upper layers of the disc, strong magnetic fields tend 
to develop during the simulation. Because of the low density 
there, the associated Alfven velocity becomes very large and 
the time step consequently very small. To prevent this from 
happening, we used the Alfven speed limiter first introduced 
by [Miller & Stone (2000) whose effect is to prevent the Alfven 
speed becoming significantly larger than a user-defined thresh- 
old v^. In the simulations presented in this paper, we used a 
uniform value of = 0.7. 

3.2. Cylindrical models 

In order to examine the effect of stratification, we ran a non 
stratified cylindrical model, labelled CI. The set-up was very 
similar to that for model S2. The vertical computational domain 
covered -0.21 < Z < 0.21, and the initial magnetic field was a 
zero-net flux toroidal field defined by: 

B^ = Bocos(67r^— ^1 . (21) 

\ «o»/ - "in / 

Bo is calculated such that the average value of < /3 >= 25. We 
also consider another cylindrical disc run, C2, which used a net 
flux toroidal field with P = 270 and was described in detail 
inlpromang & Nelson (2005). Because of the vertical periodic 
boundaries, the magnetic flux is conserved in these models. 

3.3. Model descriptions 

The detailed properties of the models we ran are described in 
table [2 Column 1 gives their label. Column 2 indicates the 
value of the parameter cq. Given that GM = 1 in our simu- 
lations, this is also equal to H/R. The resolution, the size of the 
computational box in the <p and directions are respectively 



given in columns 3, 4 and 5. The type of boundary conditions 
we use at the inner and outer radial boundary and in the direc- 
tion are shown by columns 6, 7 and 8. Finally, the code we use 
for each particular model is indicated by column 9. Note that 
runs CI and C2 are cyUndrical disc models (see section lT2t . 

4. Results 

We now present the results of our simulations. Although differ- 
ences in detail are found when comparing our various runs, a 
common picture of the early evolution of our models emerges. 
The magnetised regions of the initial discs become unstable to 
the MRI on the local dynamical time scale, in broad agreement 
with linear theory. Dynamo action associated with the MRI 
causes the local field to amplify, and as turbulence develops the 
field buoyantly rises from the regions near the midplane toward 
the disc surface. Initially the stresses associated with the field 
as it rises into these upper regions cause rapid transport of an- 
gular momentum there, allowing the magnetised fluid to spread 
rapidly through the disc in the radial direction near the disc 
surface. Contemporaneously the growth of the MRI through- 
out the magnetised core of the models causes the disc to be- 
come globally turbulent on a time scale of ~ 100 orbits. The 
end result is a disc whose global structure consists of a mag- 
netically subdominant core within ~ 2H of the disc midplane 
which remains highly turbulent due to continuing instability 
to the MRI, above and below which reside a magnetised and 
highly dynamic corona which becomes magnetically dominant 
near the disc surface and stable to the MRI. Angular momen- 
tum transport and associated mass flow through the disc allows 
magnetic field to diffuse into the resistive regions near the ra- 
dial boundaries, where the fluid remains non turbulent because 
of the resistivity and viscosity there. 

4.1. Dependance on resolution 

The picture presented above is true for all our simulations in 
their early phases. However, when designing useful and accu- 
rate stratified and turbulent protoplanetary disc models, several 
issues need to be addressed before presenting any detailed sci- 
entific results. 

The duration of the models themselves is one of them. 
In global cylindrical simulations, Papaloizou & NelsoJ ( l2003l) 
demonstrated that long integration times, over several hundred 
orbits are required for the angular momentum transport quan- 
tities to reach meaningful saturated values. This is why, in this 
paper, we ran each model for at least 450 orbits, and in some 
cases in excess of 600 orbits. In the following, we will show 
that this is sufficient to reach a steady state in the underlying 
disc structure, which is established after between 250 to 500 
orbits. 

A second issue is to determine the minimum resolution re- 
quired to allow such models to maintain turbulence over long 
run times. To do so, we present here the results of models 
Sla and Sib. They were respectively run with GLOBAL and 
NIRVANA and have a resolution (A^„ A^^, Ng) = (272, 90, 126), 
corresponding to 15 grid cells per scale height in the vertical di- 
rection (^ 8 zones per scale height in the azimuthal direction. 
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Fig. 2. The ratio of the wavelength of the most rapidly growing MRI mode A,,, to the local cell spacing A defined in the text. From 
left to right, the different panels show results from models Sla, Sib and S2. Any region coloured white has Am/ A < 5, indicating 
that the fastest growing mode of the MRI is not well resolved there. 
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Fig. 1. Time history of (um) obtained in model Sla (solid) and 
Sib (dashed line). The agreeement between the two curves 
shows that GLOBAL and NIRVANA produce similar results. 
However, in both runs, the rate of angular momentum trans- 
port displays signs that it decays with time, indicating that the 
resolution is too low in models Sla and Sib. At the end of 
both simulation, angular momentum transport is weak, with 
(um) ~ lO"-' in both models. 



^ 3 zones per scale height in the radial direction at the disc in- 
ner edge and ^ 22 zones per scale height at the disc outer edge). 
Being identical in their set-up, these two models are also useful 
as a direct comparison between the two codes. For both cases, 
we find that the MRI grows leading to fully developed MHD 
turbulence and outward angular momentum transport. 

The time history of (am) is shown in figure[2for model Sla 
(solid line) and Sib (dashed line). The two curves show very 
good agreement overall. This simply indicates that GLOBAL 
and NIRVANA give very similar results for the same prob- 
lem. It is therefore meaningful to compare the results obtained 
by the two codes when they use different starting conditions 
and physical parameters. However, the slow decrease of (qm) 




200 300 
Time (in orbits) 



Fig. 3. Time history of the toroidal magnetic flux F^(0o,t) 
threading the disc, normalized by its value at f = 0, obtained 
in model S2. The solid line corresponds to Oq - 0.3, the dotted 
line to 0() = 0.2 and the dashed line to Oq = 0.1. All of them 
show that the flux is expelled out of the midplane of the disc 
toward its corona. 



with time in both models shows that MHD turbulence is get- 
ting weaker as the simulations proceed. No steady state seems 
to be reached in these simulations, apparently because the res- 
olution used in models Sla and Sib is insufficient for vigorous 
MHD turbulence to be sustained in long runs. Indeed, in order 
to maintain turbulence, the simulations must be capable of re- 
solving the unstable modes of the MRI. The critical wavelength 
for instability is given by (Balbus & Hawley 1991): 



Ar 



27T Va 

V3 n 



and the wavelength of the fastest growing mode is 




(22) 



(23) 



S.Fromang Sl R.P.Nelson: Stratified protoplanetary discs models 



7 



where is the Alfven speed defined by 

A simulation must have at least 5 grid cells per wave- 
length Am for the fastest growing mode to be resolved (e.g. 
iHawlev et al]ll995t iMiller & Sto ne 2000). Figure El shows a 
contour plot in the (r, 6) plane which displays the ratio A„/A 
for models SI a. Sib and S2. Here A is the local cell spac- 
ing (the diagonal distance between cell vertices along a path 
that passes through the cell centre). is calculated from az- 
imuthally and temporally averaged values of the magnetic field, 
B, and density p. The time averages were performed between 
the 500* and 600* orbit of models S la and Sib. Regions of 
figurelllwhich are coloured white correspond to regions where 
Am/ A < 5. It is clear that the MRJ modes with wavelengths be- 
tween Ac and /!„, are not well resolved in these models around 
the midplane, leading to the weak and declining angular mo- 
mentum transport shown in figure^] 

Motivated by these results, we increased the resolution by 
a factor of 5/3 to 25 vertical grid cells per scale height. This 
translates into a resolution {Nr,Nii„Ng) - (455, 150,213). As 
we demonstrate in section l4~2l the results of this model, la- 
belled S2, suggest that a steady state is reached after about 300 
orbits (see for example figure|5|l for this larger resolution. The 
rate of angular momentum transport seems to saturate, giving a 
saturated value of (a) ~ 4 x 10"^ for the remainder of the sim- 
ulations. A contour plot of A,„/A for model S2 is shown in the 
right panel of figure |2 where A,„ was time averaged between 
the 350* and 450* orbits. It is clear that A,„ is well resolved 
throughout the disc in this model, leading to the turbulence be- 
ing sustained. The value of A„ in the corona exceeds the total 
disc model thickness of ^ 8.5// for values of \9-7t/2\ > 3H/R, 
such that the corona is stable against the MRI. In the rest of 
this paper we will describe results of models obtained using the 
same resolution as model S2. Unfortunately, this leads to very 
long computing times, at the limit of present day capacities: 
on standard Pentium 2.8GHz Xeon chips each model requires 
about 50000 CPU hours, or ~ 6 CPU years. Of course, it is still 
possible that MHD turbulence dies on time scales of several 
hundreds of orbits even for this large resolution, but the limited 
computational power available at the present time prevents run- 
ning them for longer than about 500 orbits. It is probably also 
the case that these simulations have not reached full numerical 
convergence. Figure|2indicates that the smaller wavelength un- 
stable modes ~ Ac are at best marginally resolved in model S2, 
such that a larger resolution may allow these modes to be more 
active. Nevertheless, the saturated state we generally obtained 
for the last 200 to 300 orbits sustains turbulence well enough 
to extract meaningful diagnostics describing the structure of the 
disc. 

It is rather surprising that model Sla and Sib fail to show 
sustained MHD turbulence, as the resolution used is similar 
to that used by Nelson (20ol) and|^omang & Nelson (2005), 
who reported sustained turbulence in cylindrical disc models 
initiated with net flux toroidal magnetic fields. This is because 
the magnetic field is trapped in the midplane of the disc in these 



cylindrical models while it is expelled to the coronae and out of 
the computational domain through the open boundaries in the 
stratified models presented in this paper. To illustrate this result, 
we plot in figure|3the time history of the toroidal magnetic flux 
as obtained in model Sla. It is defined as 

F^OoJ)^ B^(r,Q,0,t)rdrd0 , (25) 

and in figure|3lis normahzed by its initial value - F^{0o, 0). 
The different curves in this figure corresponds to Oq - 0.3 (solid 
line), 6q — 0.2 (dotted line) and 6q — 0.1 (dashed line). All de- 
crease from 1 to small values during the first 100 orbits (which 
corresponds to the growth of the MRI throughout the whole 
disc) and then oscillate around zero. This behavior shows that 
the magnetic flux is gradually expelled from the midplane to- 
ward the corona of the disc and out of the computational do- 
main (the delay shown by the solid line to decay adds further 
weight to this conclusion). The oscillations between positive 
and negative values indicate that the adopted (vertical) bound- 
ary conditions allow azimuthal net flux into the domain, al- 
though essentially zero poloidal net flux enters. In the midplane 
of the disc, the amplitude of the oscillations is about 20 - 30%. 
After about 150 orbits, the disc midplane essentially behaves 
as if MHD turbulence had been initiated using a zero net flux 
toroidal magnetic field. This requires a larger numerical reso- 
lution for the turb ulence to be s ustained over large periods (see 
the discussion bv lNelsonll2005l) and explains why models Sla 
and Sib fail to reach a saturated state while cylindrical models 
using an equivalent resolution do. 

4.2. The fiducial run - model S2 

We describe in this section the general properties of model S2, 
which was run for 500 orbits. As described above, the MRI 
grows in the first 100 orbits, before developing into MHD tur- 
bulence, which then diffuses over the entire computational box 
(except for the buffer zones described in section im . The disc 
settles into a quasi-steady state after 250 orbits for the remain- 
der of the simulation. The structure of the disc after 400 or- 
bits is illustrated in figure|3 The left hand snapshot shows the 
distribution of the logarithm of the density in the (r, 6) plane 
and the right hand panel represents the logarithm of the Alfven 
speed va. Strong radially propagating density waves are seen in 
the former, while the small scale structure of the Alfven speed 
in the latter confirms the turbulent nature of the flow, and also 
shows that the disc forms a structure consisting of a turbulent, 
magnetically subdominant 'core' near the midplane where the 
Alfven speed is small, and a magnetically dominant corona in 
the upper regions of the disc where the Alfven speed is large. 

4.2.1 . Angular momentum transport 

In this section, we quantify the rate of angular momentum 
transport resulting from the MHD turbulence. It is measured 
by the sum of the Maxwell and Reynolds stresses, akeady dis- 
cussed in section ITTI The time history of (a), (om) and {ur) 
are plotted in figure Irrespectively with the solid, dashed and 
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Fig. 4. The left panel shows the logarithm of the density distribution in the (r, 0) plane for model S2. Before taking the logarithm 
we transformed p using the following: p ^ px r'-^ x exp [z^/(2 x 0.07^)j x exp [-Z2/(2x0.l2)]. This was done to increase the 
contrast in the figure. The right panel shows the logarithm of the Alfven speed in the (r, 6) plane, obtained in model S2 after 400 
orbits. 
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Fig. 5. Time history of the volume averaged value of {a) (solid 
line), (om) (dashed line) and (a^) (dotted line) in model S2. 
After an initial evolution driven by the linear growth of the 
MRI between and 250 orbits, the disc settles into a quasi 
steady state for the remaining part of the simulation during 
which (a) ~ 0.004. 



dotted l ines. As repor t ed before for global MHD disc simu- 
lations iH awlev' 2000, 20ol LSteinacker & Papaloizoul l2002t 
IPapaloizou & Nelson. .20031) . the Maxwell stress dominates 
over the Reynolds stress in the approximate ratio of 3:1. The 
rate of angular momentum transport saturates after 250 or- 
bits and shows an almost constant value of (a) ~ 4 x lO"-' 
for the remainder of the simulation. This is very similar to 
results reported for zero-net flux numerical simulations of 
cylindrical discs which further support the conclusion that 
our stratified disc model behaves Hke a zero-net flux disc. 




radius 



Fig. 6. Radial profile of a (solid line), om (dashed line) and 
Or (dotted line) in model S2. The curves are time-averaged 
between f = 350 and f = 450 orbits. 



IPapaloizou & NelsonI (l2003h report values of (a) in the range 
0.002 to 0.005 for a suite of zero net flux cylindrical disc 
models, and model CI achieves a saturated state with (or) ^ 
0.002. This slightly smaller value of (a) is consistent with the 
smaller vertical resoluti on used in model CI compared with 
the models presented bv IPapaloizou & NelsonI ll2bo.3.) . The ra- 
dial profile of the Maxwell and Reynolds stresses, normal- 
ized by P, are represented in figure |6l respectively with the 
dashed and dotted line. The solid line simply represent s a, th e 
sum of the two. As reported by Papaloizou & Nelsor| ll2003h . 
we found large fluctuations in snapshots of these quantities. 
The smooth profile shown in figure |6l results from an averag- 
ing between 350 and 450 orbits. We note that time averaged 
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Fig. 7. Vertical profile of the Maxwell stress {solid line) and 
the Reynolds stress {dashed line) obtained in model S2. The 
results are time averaged between 350 and 450 orbits and in 
radius between r - A and r - 5. Both curves are normalized by 
the midplane pressure. A decrease of both stresses is apparent 
forZ/// > 3. 



a profiles in global MHD simulations (including these) often 
show va riations between a factor of 2 and 3 across the disc 
(see e .g. IPaoaloizou & Nelsonll2003t ISteinacker & Paoaloizoul 
l2002h . The origin of these variations is not clear, but may be re- 
lated to the fact that fluctuations in a anticorrelate with changes 
to the local density and pressure due to mass transport. Local 
variations in a may thus be reinforced. It is possible that these 
variations will decrease with longer run times as the system 
loses memory of its initial conditions. 

4.2.2. Structure of the corona 

The results described above regarding the volume integrated 
properties of angular momentum transport show little differ- 
ence compared with zero-net flux cylindrical disc models. In 
this section, we detail the structure of the upper layers of the 
disc. 

Figure shows the vertical profile of the Maxwell {solid 
line) and Reynolds {dashed line) stresses, normalised by the 
midplane pressure. Both curves are averaged in space between 
r - A and r - 5 and in time between t - 350 and t - 450 orbits. 
In agreement with the results obtained in local disc simulations 
jMiller & Stondl2000l) . both are fairly constant for |Z| < 2.5H 
before decreasing in the upper layers of the disc. This change 
is due to the establishment of a strongly magnetised corona. 
This is illustrated by figure |8] which shows the vertical pro- 
file of the ratio Pmag/P at r - 3.5 (the solid line corresponds 
to model S2). Simulation data were averaged in time between 
350 and 450 orbits to smooth out the turbulent fluctuations (the 
other curves represented in figure |S] plot the results of models 
S3, S4 and S5 and will be discussed later). The relative strength 
of the magnetic field increases with Z, approaching equiparti- 
tion in the neighborhood of the lower and upper boundaries. 



100.000 F 



1 0.000 r 




0.001 I ... I ... I .. . 1^":' . . I . . . I . . , 

-6 -4 -2 2 4 6 
z/H 

Fig. 8. Vertical profile of the radio P„,ag/P at r = 3.5 obtained 
in model S2 {solid line, averaged between 350 and 450 orbits), 
model S3 {dashed line, averaged between 270 and 310 orbits), 
model S4 {dotted line, averaged between 500 and 600) and 
model S5 {dotted-dashed line, averaged between 500 and 600 
orbits). For all models, the disc is composed of a weakly mag- 
netised core and a corona where the field is close to equiparti- 
tion. The differences between the different curves results from 
the boundary conditions and are further discussed in the text. 
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Fig. 9. The time averaged values of \B,\/B, \Bg\/B and \B^\/B as 
a function of 9/{H/R) for model S5. This was obtained by time 
averaging for 100 orbits at radius r = 3.5. The (p component 
of B is represented by the solid line, the r component by the 
dotted line, and the component by the dashed line. 



Even though the flow in the corona is not unstable to the MRI, 
we found it to be highly dynamic, exhibiting strongly time- 
dependent behaviour. Nevertheless, the structure of the corona 
is quite different to that in the midplane, with much smaller 
scale magnetic field fluctuations near the equatorial plane than 
in the corona of the disc. To try to quantify the topology of the 
field in the disc, figure|9]shows the variation of \Bi-\/B, \B^\/B, 
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Fig. 10. Radial profile of the velocity fluctuations obtained in 
model S2. The curves, time-averaged between t - 350 and t - 
450 orbits, correspond to radial {solid line), azimutal (dashed 
line) and vertical {dotted line) velocity fluctuations. 



Fig. 11. Vertical profile of the velocity fluctuations obtained 
in model S2 at r = 4.5. The curves, time-averaged between 
t - 350 and f = 450 orbits, correspond to radial {solid line), 
azimutal {dashed line) and vertical {dotted line) velocity fluc- 
tuations. 




and \Bg\IB versus ZjH, where B = ^B^ + + B^. This plot }, 

corresponds to a radial location r = 3.5 and was obtained by 

time averaging for 100 orbits between t - 350 and 450 orbits. i 

It shows that the magnetic field topology is dominated by the 

azimuthal component of the field, as expected because of the ^ 

strong Keplerian shear. However, as one moves away from the 

midplane toward the disc surface, there are some indications 

that the r and 9 component of the field become more important "■ ° 

as their relative magnitude increases by about a factor of two. 

We comment here that the boundary conditions imposed on the -i 

magnetic field at the upper and lower disc surfaces were zero 

gradient outflow conditions. 



4.2.3. Velocity and density fluctuations 

It is of interest to look at the velocity and density fluctuations 
generated by the turbulence in these models. Within the context 
of pl anetary formation, they affect the evolution of dust parti- 
cles dpromang & Papaloizoui.200 6). whose spatial distribution 
is important for the observational properties of protoplanetary 
discs, as well as the orbits of lar ger bodies such as boulders, 
planetesimals and protopla nets (iNelspn & Papaloizoui i2004t 
lNelsonl2005t iRom ang & N elsonl2005h . 

Figure shows the radial profile of the velocity pertur- 
bations obtained in model S2, averaged between t - 350 and 
t - 450 orbits. The solid, dashed and dotted line respectively 
represents the radial, azimuthal and vertical velocity pertur- 
bation, normalized by the local sound speed and averaged in 
space within ///2 around the disc midplane. Typically, these 
fluctuations are all of the order of a few percents of the sound 
speed. This is sligh tly smaller than previo us results obtained 
in the shearing box dStone et al.lll996t iFro mang & PapaloizoiJ 
I2006 '). who found values of the order of 10%. There is a marked 
tendency for the radial fluctuations to be larger, by about a fac- 
tor of two, than the azimutal and vertical fluctuations. This is 



a i i ^ " 

Fig. 12. Snapshots of the Mach number M, in the {r, 0) plane at 
t - 450 orbits obtained in model S2. Motions in the bulk of the 
disc are largely subsonic while weak shocks (with Ms between 
1 and 1 .5) can be seen in its corona. 



due to the presence of density waves travelling radially through 
the disc. This is also seen in local boxes but its efi'ect is less 
pronounced in that case. The vertical profiles of the fluctua- 
tions are shown, using the same con ventions, in figurelll| As 
noted by lMiller & Stonel tOOdi and lTurner et~lj (.2006) who 
presented shearing box calculations, they increase in the upper 
layers of the disc, where the averaged Mach number reaches 
0.4. This increase in perturbed velocities arises because distur- 
bances that are excited near the high density midplane increase 
in amplitude as they propagate vertically into the lower density 
regions, due to conservation of wave action. In addition, the 
increasing influence of the magnetic field with height means 
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Fig. 13. Normalised frequency distribution of the relative den- 
sity fluctuations in models S2 {solid line), S5 (dashed line), CI 
(dot-dashed line) and C2 (dotted line). All curves are averaged 
in time for 100 orbits. Results obtained in run S2 and S5 are av- 
eraged in the volume |Z| < H and 3 < r < 5. Results obtained 
for runs CI and C2 are averaged over the entire vertical extend 
and in radius between r = 3 and r = 5 for model C 1 , and be- 
tween r = 2.5 and r - 4 for model C2. The curves for runs 
S2, S5 and CI have a similar width, while model C2 produces 
a wider distribution because it contained a net flux magnetic 
field. 



that propogating Alfven waves in the disc corona can excite ap- 
proximately sonic motions where the Alfven speed exceeds the 
sound speed. The result is that shocks are generated in the up- 
per regions of the disc, as illustrated by figure [Qlwhich shows 
a snapshot in the (r, 6) plane of the perturbed velocity divided 
by the sound speed. Weak shocks, for which M, reaches 1.5, 
are visible on this image. 

The normalised frequency distribution for the density fluc- 
tuations is shown in figure [O] for model S2 (solid line). The 
results of the simulations were averaged in time between 350 
and 450 orbits in order to smooth out the turbulent fluctuations. 
An additional spatial averaging was performed in the volume 
3 < r < 5 and |Z| < H. The results show that the distribution 
is approximately Gaussian with standard deviation a-p - 0.08. 
This result agrees well with those of models S5 (dashed line) 
and CI (dot-dashed line), showing that neither the boundary 
conditions nor the stratification have a strong impact on the am- 
plitude of the turbulent density fluctuations near the midplane. 
However, the results of model C2, shown with the dotted line, 
produce a broader distribution (a-p ^ 0.13). This is because the 
toroidal net flux is nonzero in this model, while it vanishes in 
model S2, S5 and CI. This results in a stronger MHD turbu- 
lence, along with slightly larger density fluctuations. It is al- 
ready known that these density fluctuations may be responsible 
for driving sto chastic migration of low mass planets and plan- 
etesimals (e.g. iNelson & Pa oaloizou 2004: Nelson 2005), and 
this effect will be explored in more detail in a forthcoming pa- 
per within the context of stratified, turbulent disc models. 
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Fig. 14. Time history of the Maxwell stress in model S2 (solid 
line) and S3 (solid line). The latter has a larger vertical extent. 
The good agreement between the two curves shows that the 
boundary conditions and vertical extent have little effect on the 
results of model S2. 



4.3. The effect of the boundary conditions in 9 

We now present the results of simulations designed to test the 
role of the vertical boundary conditions we have adopted. We 
begin by examining a model whose vertical domain is larger 
than model S2 but otherwise similar, and then describe a model 
with periodic boundary conditions in the vertical direction. 

4.3.1 . A larger vertical domain - model S3 

We first focus on model S3. Compared to model S2 described 
above, it has an extended domain in 0: \0 - 7r/2| < 0.4. This 
corresponds to a total of 5.7 scale heights on both side of the 
equatorial plane. For the actual resolution to remain the same, 
the total number of grid cells was increased to (Nr,N^,Ng) = 
(455, 150, 285). AH the other parameters of model S3 are ofli- 
erwise identical to those of model S2. Note however that the 
larger number of cells increased the computing time of model 
S3. As a result it was run for only 325 orbits. 

Figure compares the time history of (qm) in model S2 
(solid line) and in model S3 (dashed line). Both curves show 
similar evolution. Moreover, model S3 starts to saturate after 
about 270 orbits at a level that is comparable to that of model 
S2. Figure 1^1 also shows the vertical profile of P„mglP (dotted 
line), averaged between 270 and 310 orbits. It shows that model 
S2 and S3 agree very well in the bulk of the disc. 

These results demonstrate that the structure of the discs de- 
scribed so far are not affected strongly by the boundary condi- 
tions and vertical extent of the computational domain. 

4.3.2. Periodic boundary conditions - model S4 

We have computed one model, run S4, which uses periodic 
boundary conditions in the vertical direction (i.e. 0-direction). 
Clearly this is not a realistic boundary condition to use for an 
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Fig. 15. Time evolution of the volume averaged stress parame- 
ter (qtm) for model S4 (solid line) and S2 (dashed line). Despite 
the different boundary conditions, the two models display simi- 
lar time history for the angular momentum transport properties. 
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Fig. 16. The time evolution of the azimuthal magnetic flux for 
model S4, normalised to its initial value. The dashed line shows 
the total flux within an angular distance - n/2\ < H/R above 
and below the midplane. The dotted line shows the flux within 
an angular distance |0-7r/2| < 2///7? above and below the mid- 
plane. The solid line shows the flux throughout the whole disc, 
and varies with time due to the imperfect nature of the periodic 
boundary conditions adopted in the meridional direction. 



accretion disc, but it nonetheless provides a test of the role that 
the vertical boundaries play in these simulations. Before de- 
scribing the results of the simulation in detail, we first comment 
that the periodicity condition implemented in the NIRVANA 
(or GLOBAL) code is imperfect when applied in the merid- 
ional direction. This is because the physical sizes of the cells 
that overlap at the top and bottom of the disc are different. The 
effect of this is small, but one manifestation is that magnetic 
flux is not conserved as it passes through the upper surface 
of the disc and re-enters through the lower surface (and vice 
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Fig. 17. A slice plotted in the (r, 6) plane of the disc showing 
the logarithm of the Alfven speed for model S4. Comparing 
with similar plots for models S2 (figure |3 and S5 (figure |22j 
shows that the trapping of magnetic flux in the upper disc re- 
gions of model S4 leads to an elevated value of the Alfven 
speed there. 



versa). As shown below, this causes a 5 % decrease in mag- 
netic flux in the domain over the simulation run time. 

The time evolution of the volume averaged (aw) value is 
shown in figure for model S4 (solid line) and compared to 
model S2 (dashed line). The early rise and fall of this quantity 
is similar in both models. One feature of the runs that appears 
throughout this paper is that the time required for saturation of 
the turbulence is longer for the NIRVANA simulations than for 
the GLOBAL ones. Nonetheless the final states that are reached 
are very similar in each case. Inspection of figure [Tsl shows 
that during the time interval 350-550 orbits, short duration out- 
bursts are observed in the um values. This phenomenon seems 
to be explained by figure ^] which shows the time evolution 
of the azimuthal magnetic flux in the computational domain de- 
fined by equation ( 125 > . The solid line shows the total flux in the 
domain, and because of the periodicity of the vertical bound- 
aries this is approximately conserved (the small non conserva- 
tion was explained above). The dashed and dotted lines show 
the magnetic flux in the disc below heights \6 - 7r/2| < H/R 
and \6 - 7t/2\ < 2H/R, respectively. As described for model 
S2, the onset of the MRI and turbulence causes the magnetic 
flux to rise up through the disc to form a magnetically dom- 
inated corona. In model S4 this magnetic flux is effectively 
trappped in the disc and is able to build up to large values in 
the corona. Figure[T6lshows that substantial flux occassionally 
comes down from the corona into the turbulent core, within 
~ 2H/R of the midplane, where its presence causes an episodic 
increase in turbulent activity. 

The effects of trapping the magnetic flux in the domain on 
the evolution of the disc is illustrated by figure[n]which shows 
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Fig. 18. Snapshot of velocity fluctuations in the disc normalised 
to the local sound speed, plotted in the (r, 6) plane. It is clear 
that strong shocks (Mach numbers > 10 ) are generated in the 
upper disc regions due to the large magnetic forces there. 



a snapshot of the Alfven speed plotted as a slice in the (r, 6) 
plane. The contrast in the Alfven speed between the inner disc 
core and the upper corona is much greater in this case than was 
observed for model S2 because the magnetic field in the upper 
regions is larger. This is further illustrated by figure |8] which 
shows the ratio of magnetic to thermal pressure Pmag/P as a 
function of disc height for model S4 using the dotted line. This 
plot corresponds to r = 3.5 and was obtained by time averaging 
between t - 500 and t - 600 orbits. It is clear that the corona 
in this case is highly dominated by the magnetic field, with 
Pmag/P - 30 near the disc surface, in contrast to the value of 
Pmag/P - 1 obtained in model S2. 

One eff'ect of these large magnetic field and Alfven speed 
values is that large supersonic motions are generated in the disc 
corona. This is illustrated by figure[^which shows a snapshot 
of the velocity perturbation divided by the local sound speed, 
projected onto the (r, 6) plane. Strong shocks can be seen in 
the upper disc surface with Mach numbers M, in excess of 10 
being common. This is in contrast to the more quiescent state 
of the corona obtained in model S2 where Mach numbers be- 
tween 1-2 are observed. We note that M, ~ 10 indicates that 
the maximum gas velocities in the corona are basically given 
by the Alfven speed limiter (which is roughly ten times the 
sound speed). This is probably a numerical consequence of the 
particular value chosen for y^. However, the exact strength of 
these shocks is not of crucial importance since they are mostly 
a consequence of the unphysical periodic boundary conditions, 
which serve to illustrate the effect of trapping the toroidal net 
flux on the structure of the corona. 
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Fig. 19. Time history of the stress parameter (a) for model S5. 
The dotted line shows the contribution from the Maxwell stress, 
the dashed line shows the contribution from the Reynolds 
stress, and the solid line shows the sum of these. 



4.4. A thicker disc - model S5 

We now present the results of run S5 whose parameters are de- 
scribed in table[2 This run used a thicker disc with H/R = 0.1, 
and a correspondingly smaller number of grid cells (A^^, A^^, 
A^e)=(360, 120, 210). A distinct advantage of using a thicker 
disc is that in principle it allows the use of a smaller number 
of grid cells while still giving rise to a resolved model. We re- 
mind the reader that model S5 used a "viscous outflow bound- 
ary condition" at the inner edge of the computational domain 
(see section lTTI for details). 

4.4.1 . Global disc properties 

The early evolution of this model proceeds very much along 
the lines discussed in the opening paragraph of section|3] The 
time evolution of the volume averaged viscous stress parame- 
ter (a) is presented in figure[^ and shows that it saturates at a 
value of (a) ^ 4x lO"-' after 500 orbits, in basic agreement with 
the models S2, S3, and S4. The time evolution of the azimuthal 
magnetic flux defined by equation ( I25t is shown in figure |20| 
Similar behaviour is seen in model S5 as was described for run 
S2, with magnetic buoyancy causing the initial flux to rise ver- 
tically through the disc into the corona where it escapes through 
the open boundary. Close inspection of the dot-dashed line in 
figure|20|shows that the disc near the midplane contains essen- 
tially zero net magnetic flux after about 100 orbits, such that 
sustained MHD turbulence there requires the action of a dy- 
namo. As already described, this feature of these simulations 
adds to the computational burden as zero net flux simulations 
require high resolution in order that numerical resistivity does 
not quench the MRl. The solid line in figure|20|shows the total 
flux in the whole computational domain. The fact that it os- 
cillates about the zero line indicates that the adopted boundary 
conditions lead to magnetic flux entering the computational do- 
main. The relatively small amount of flux that enters suggests 
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Fig. 20. Time history of the azimuthal magnetic flux for model 
S5. The dot-dashed line shows the normalised magnetic flux 
contained within \9 - 7r/2| < HjR about the midplane, the 
dashed line shows the flux within \9 — 7r/2| < 2H/R, the dot- 
ted line shows the flux within \6 - 7r/2| < 3H/R and the solid 
line corresponds to the magnetic flux within the whole disc. As 
for model S2, the initial flux within the disc escapes through 
the upper and lower disc surfaces leading to a substantial re- 
duction throughout the disc. In particular, the regions around 
the midplane tend towards a state of zero net flux. Oscillations 
in the total magnetic flux are caused by magnetic flux entering 
through the vertical boundaries. 



that this feature does not dramatically alter the results of our 
simulations. 

Figure]^ shows snapshots of the density in the disc after 
500 orbits. The left panel shows the density in the disc mid- 
plane, and the right panel shows a vertical slice through the (r, 
ff) plane. The radial transport of mass during the simulation has 
caused a shallow depression to form in the density profiles at 
radii between 5 < r < 7 which are apparent in the figure. Also 
apparent are the trailing spiral waves excited by the turbulence. 
These propagate radially through the model and, because the 
disc is isothermal in the vertical direction, these waves propa- 
gate with little vertical structure. Animations of the disc density 
projected onto the (r, 6) plane indicate that individual promi- 
nent spiral waves which propagate radially occupy most of the 
vertical extent of the disc. Comparison between figures 1211 and 
13 shows that the 'viscous outflow' boundary condition at the 
inner edge of the disc is having the desired effect of preventing 
a substantial build-up of mass there. As we comment in sec- 
tion l4.4.3l some build-up of mass near the inner boundary does 
occur in this model because the chosen value of the outflow ve- 
locity was too small, but this should be a simple problem to 
remedy in future models. 

Figure |22lpi"esents a vertical slice projected onto the (r, 6) 
plane showing the logarithm of the Alfven speed. As observed 
in the runs S2 and S4, the Alfven speed noticeably increases as 
one moves away from the midplane to the disc surface above 
a height of about 2H. It is also clear from this figure that the 
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Fig. 22. This figure shows the logarithm of the Alfven speed in 
the disc after 500 orbits obtained by plotting a slice in the (r, 6*) 
plane. 



disc is magnetically active all the way down to the inner radial 
boundary due to the 'viscous outflow' condition used in this 
model. The increase in relative strength of the magnetic field 
with height is also illustrated by figure |S] where the ratio of 
magnetic pressure to gas pressure is plotted as a function of 
at radius r - 3.5 using the dot-dashed line. In agreement with 
the results of model S2, one sees that Pmag/P - 0.01 within 
\d - n/2\ < 2H/R, but it quickly rises to P,„ag/P ^ 0.5 in the 
low density regions above the midplane. 

The variation of stress with height is shown in figure 
which shows a similar fall-off" in the Maxwell and Reynolds 
stresses with height as observed for model S2. This figure cor- 
responds to radius r = 3.5, and was obtained by time averaging 
the stresses for 100 orbits. 

Figure|23 which should be compared with figure|9l shows 
the variation of \Br\/B, \B^\/B, and |Bfl|/B versus 6. This plot 
corresponds to radial location r = 3.5 and was obtained by time 
averaging for 100 orbits between f = 500 - 600 orbits. It shows 
a bigger increase of the r and 6 component of the field in the 
upper layers of the disc than in model S2. This is clearly influ- 
enced by the boundary conditions imposed at the disc surface, 
where the field is defined to be normal to the boundary with 
magnitud e defined by the V.B - condition in the NIRVANA 
runs (e.g. lHawlevl l2000). but, as suggested by the results of 
model S2 (which show a similar but reduced effect with differ- 
ent boundary conditions), it is probably also due to the vertical 
stretching of field lines as localised regions of magnetised fluid 
rise up from the disc midplane. 

4.4.2. Velocity and density fluctuations 

The ratios of the velocity fluctuations in the r, 6, and <p direc- 
tions to the local sound speed were calculated and found to 
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Fig. 21. The left panel shows an image of the density at the disc midplane after 500 orbits for model S5. The right hand panel is a 
slice plotted in the (r, Q) plane showing the vertical density profile (this panel shows the logarithm of the density after performing 
a transformation of p similar to that described in the caption of figurel^. These images show that the action of the stresses have 
caused signification mass transport within the disc, resulting in an apparent density depression at r ~ 6. 



be in very good agreement with the results of run S2 plotted 
in figure ^2 The increase of the velocity perturbation Mach 
number as one moves from the midplane to the disc surface 
results in weak shocks being generated in the corona, as illus- 
trated by figure l25l which shows a snapshot in the (r, 6) plane 
of the perturbed velocity divided by the sound speed. As ob- 
served in model S2, and unlike model S4, typical Mach num- 
bers for these shocks range between 1-3. The fact that model 
S4 showed much stronger shocks illustrates the role that mag- 
netic forces in the corona have in exciting these supersonic mo- 
tions. 

The distribution of the density fluctuations in model S5 is 
represented in figure [21 with the dashed line. Simulation data 
were averaged in time between 500 and 600 orbits, and in space 
in the radial interval 3 < r < 5 and in the height interval \0 - 
7r/2| < HjR to obtain this curve. Once again, the results are in 
very good agreement with those of model S2. 

4.4.3. Mass flow in the disc 

Once model S5 had completed just over 600 orbits it was 
stopped. It was restarted again at the 500 orbit mark, but with its 
density field and azimuthal velocity reset to the values they had 
initially at time t - 0. All other variables (e.g. magnetic field, 
vertical and radial velocities) had the values corresponding to 
the 500 orbit mark of the S5 run. This procedure was under- 
taken to smooth out the large variations in surface density that 
arise during the early stages of these simulations because the 
stresses have large radial and temporal variations during these 
times. The aim is to generate a turbulent protoplanetary disc 
with simple surface density structure as a function of radius. 

Having restarted the model it was run for 55 orbits, at which 
point time averaging of the state variables and stresses within 
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Fig. 23. The time averaged values of the Maxwell stress {solid 
line) and Reynolds stress {dotted line), normalised to the mid- 
plane pressure, as a function of 6 for model S5. This was ob- 
tained by time averaging for 100 orbits at radius r - 3.5. This 
figure is similar to that obtained for model S2, showing a simi- 
lar drop-off" of stress with height. 



the disc were commenced. The simulation was continued for a 
further 100 orbits while the time averaging was performed. 

The radial variation of the time averaged a values in the 
disc are shown in figure|26l The results are similar to those ob- 
tained in model S2. We note here, however, that the value of 
a is small (~ 10"^) near the disc inner radial boundary where 
a 'viscous outflow' condition is imposed on the radial veloc- 
ity. This arises because during the initial phases of the run for 
model S5, the magnetic field was non zero only in the range 
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Fig. 24. Same as figure|9] As in model S2, the azimuthal com- 
ponent remains dominant throughout the vertical domain, but 
the field topology changes near the disc surface where the 6 and 
r components become larger. 
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Fig. 26. Radial variation of the time averaged value of a for 
model S5. The curves were obtained by time averaging be- 
tween 500-600 orbits. The dashed line corresponds to the 
Maxwell stress, the dotted line to the Reynolds stress, and the 
solid line to the sum of these. 
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Fig. 25. Snapshots of the velocity fluctuations normalised by 
the local sound speed obtained in model S5. The results agrees 
very well with those of model S2 (see figureO 



2.5 < r < 6. As turbulence develops, mass and magnetic field 
are transported inward, but at a rate which is larger than ac- 
counted for by the 'viscous outflow' condition. The density in 
the midplane thus increases in this region, such that the unsta- 
ble modes of the MRl remain small and unresolved here [see 
equation ( I22H . Consequently turbulence in the inner most re- 
gion of the disc remains weak. A remedy for this would be to 
modify the 'viscous outflow' condition so that it responds more 
accurately to the inflow of matter from further out rather than 
using a prescribed inward velocity as was done for model S5. 

Figurel^shows the time averaged radial mass flux through 
the disc as a function of radius. Each line corresponds to mass 



flow at a different height in the disc. The dotted line corre- 
sponds to the region within A6 < H/R about the midplane, 
where A6 = \6 - n/2\. The dashed line corresponds to the re- 
gion bounded hy H/R < AO < 2H/R, the dot-dashed corre- 
sponds to 2H/R < AG < 3H/R, and the dot-dot-dot-dashed 
line corresponds to 3H/R < AO < (Omax or 0,„i„). The total ra- 
dial mass flux is shown by the solid line (and is the sum of 
all the other lines). Evidently negligible mass is transported in 
the upper most parts of the disc corona, as there is very lit- 
tle mass there. The other three regions considered, however, 
all contribute significantly to the mass flux. In the outer re- 
gions of the disc the zone near the midplane appears to be 
transporting mass inward whereas the upper regions of the disc 
are transporting it outward, suggesting that the long term mass 
flow in vertically stratified turbulent discs can be a complicated 
function of dis c height. A s imilar picture was described by 
lOe Villiers & HawlevI ll2003l) for simulations of accretion tori 
around black holes. 

The total vertical mass flux through both boundaries at the 
disc surface is shown in figure|28l It is clear that the disc drives 
a vertical mass flow at a rate that is less than two orders of mag- 
nitude below that which occurs radially in the disc. A similar 
result was obtained in model S2. However, the restricted size of 
our meridional domain prevents us from commenting in detail 
about any wind that may be launched from the disc surface. 

We now turn to the question of how well the averaged ra- 
dial velocity in the disc agrees with the expectations of thin 
disc theory as described by equation il5\ . We computed the 
time average of X(/?), Tm{R), Ti{{R) and using a simple finite 
difference approximation calculated the expected radial veloc- 
ity profile vs in the disc. The results are shown using the dotted 
line in figure j29L where the actual value of obtained in the 
simulation is shown using the solid line. Although the effect of 
fluctuations remain, the agreement between the predicted and 
actual values is remarkably good. This demonstrates that the 
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Fig. 27. Radial mass flux versus radius for model S5. The dot- 
ted line corresponds to \9 - 7r/2| < HjR, the dashed line corre- 
sponds to HjR <\0 - 7t/2\ < 2H/R, the dot-dashedline corre- 
sponds to 2////? < \0-7t/2\ < 3H/R, andthe dot-dot-dot-dashed 
line corresponds to 3H/R < \6 - 7t/2\ < disc surface. The solid 
line represents the total radial mass flux through the disc. 



vertically stratified turbulent discs considered here behave very 
much like standard a discs when their evolution is considered 
over long time scales. On short time scales, however, the differ- 
ences are self-evident. We note that a comparison of the pre- 
dicted and actual values of vr was undertaken for model S4 and 
gave rise to a very similar level of agreement. 

5. Conclusions 

In this paper we have presented the results of 3-D MHD sim- 
ulations of stratified and turbulent protoplanetary disc models. 
Our primary motivation is to develop disc models that can be 
used to examine outstanding issues in planet formation such as 
the migration of protoplanets, the growth and settling of dust 
grains, gap formation and gas accretion by giant planets, and 
the evolution and influence of dead-zones. Given that these 
phenomena occur on secular time scales, a key requirement is 
the development of disc models which are able to achieve a sta- 
tistical steady state and sustain turbulence over long run times. 

We examined the issue of numerical resolution, and found 
that disc models with ^ 15 vertical zones per scale height in 
the vertical direction showed a continuing slow decline in their 
magnetic activity, and gave rise to relatively small values of 
{a). A suite of higher resolution runs with ^ 25 zones per 
scale height achieved statistical steady states with values of 
(a) ^ 4 X 10""', and it was shown that these models resolve 
the fastest growing modes of the MRJ throughout the disc once 
a turbulent steady state has been achieved. For this reason we 
focused on simulations performed using this higher resolution. 
The key features of the resulting disc models are: 

- Any toroidal magnetic flux that is initially present within 
the disc is quickly expelled from the midplane due to mag- 
netic buoyancy. This occurs on a time scale of ~ 100 or- 
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Fig. 28. The sum of the mass fluxes through the vertical bound- 
ary of the disc located at9 - Omax and 9mm- 



bits, which is the time required for the MRl to grow and 
develop into non linear turbulence throughout the disc. The 
disc then evolves as if it is threaded by an approximately 
zero net flux magnetic field, such that high resolution is re- 
quired to maintain turbulent activity. 

- A quasi-steady state turbulent disc is obtained after run 
times of between 250 - 500 orbits, depending on the model. 
The volume averaged value of the effective viscous stress 
paramater {a) ^ 4 x 10"^, and time averaged radial profiles 
of a yield variations of no more than a factor of two within 
the active domains of the disc models. These results are in 
basic agreement with previous studies of cylindrical discs 
(ISawley 2001; Papaloizou & Nelson 2003) 

- The discs can be described as having a two-phase global 
structure as a function of height: a dense, magnetically sub- 
dominant turbulent core that is unstable to the MRl in re- 
gions within |Z| < 2.5H of the midplane, above and below 
which exists a highly dynamic and magnetically dominated 
corona which supports weak shocks and is stable against 
the MRl. The engine that drives this structure is the MRl 
which generates and amplifies magnetic field near the mid- 
plane, which then buoyantly rises up into the low density 
corona where it dissipates and flows out through the bound- 
aries at the disc surface. This is in basic agreement with 
the sh earing box simulations presented by Miller & S tond 
teOOd) and recent stu dies of thick tor i orbiting around black 
holes tHawlev 2000; Hawlev & Krolik.200klHawlev et alJ 
l200l[ lDe Vifliers et al. 2003). 

- The velocity and density fluctuations generated by the mod- 
els were found to be smaller than those obtained in cyUndri- 
cal disc simulati ons using tor o idal net flux magnetic fiel d 
configurations (Nelson 2 005t iFromang & Nelsoi] l2005l) . 
<5p/po - 0.08 in the stratified runs whereas in the cylindri- 
cal disc runs with net flux 5p/po - 0.13. This has implica- 
tions for the dynamics of dust, planetesimals and low mass 
protoplanets in turbulent discs as their stochastic evolution 
is driven by these fluctuating quantities. 
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Fig. 29. Comparison between the time averaged radial velocity 
obtained in the simulations (solid line) and the predicted value 
of vr obtained from equation il5i shown by the dotted line. 



- The vertically, azimuthally and time averaged values of 
the radial velocity in some of the disc models were com- 
pared with the expectations of viscous disc theory, and were 
found to give very good agreement. We conclude that, sub- 
ject to a suitable time average, global evolution of these 
stratified turbulent models is in good accord with stan- 
dard viscous disc theory (e.^. [ ^ hakura & Sunva ev , 1973 : 
iBalbus & PanaloizoiJl999HPanaloi70u & Nelsonl2003l) 

There are a number of outstanding issue s raised by our sim- 
ulation results. iFromang & NelsonI {2005) reported the pres- 
ence of anticyclonic vortices in turbulent unstratified cylin- 
drical disc models. In the stratified models we present in this 
paper, however, we did not find any evidence of vortices. 
This could be due to a number of reasons. First, the verti- 
cal stratification may pre vent the formation of vorti ces near 
the midplane. A study by Barra nco & MarcusI ( l2005l) showed 
that column-vortices in stratified discs are unstable and are 
quickly destroyed. They showed that vortices could form in 
the more strongly stratified upper regions of the disc which are 
more than one scale height above the midplane. In our simula- 
tions these regions are typically dominated by magnetic fields, 
whose associated stresses may prevent the formation of vor- 
tices there. Second, we adopted a smaller azimuthal domain 
than lFro mang & Nelson (2005): n/A- versus 7r/2 and In mod- 
els. This may prevent the formati on of vortic es, which were 
found to be quite extended in <p bv lFromang & Nelson C2005 ). 
We did not observe any vortices in the cylindrical disc model 
CI described in section l3T2l and this may be partly explained 
by the smaller azimuthal domain. Finally, the different mag- 
netic field topology conta ined in the disc may play a role: 
iFromang & NelsonI ( l2005l) used a net flux toroidal magnetic 
field. In the stratified models, the toroidal flux is quickly ex- 
pelled from the disc midplane and the evolution is more sim- 
ilar to that of a zero net flux disc. The cylindrical model CI 
contained a zero net flux magnetic field. The properties of the 
field can affect the turbulence and hence the formation of vor- 



tices. In particular, stronger spatial and temporal variations in 
the stresses may cause the surface density variations to differ 
systematically between the models presented here and those in 
f romang & Nelson (2005). If the formation of vortices in the 
models described in IFromang & NelsonI (l2005l) are related to 
the 'planet modes' described bv lHawlevl(ll987|) . then these dif- 
ferences may explain the lack of vortices seen in the stratified 
models and model CI. These issues, and their influence on the 
evolution of solid bodies will be explored in greater detail in a 
future paper. 

Finally, this is the first paper in a series which describes 
an approach to setting up models of turbulent, stratified proto- 
planetary discs capable of sustaining turbulence over long run 
times. Future papers will present a systematic study of out- 
standing problems in planet formation theory, such as disc- 
planet interactions, dust and planetesimal dynamics, and effects 
related to the presence of a dead zone. We also note that these 
models themselves can be further improved by including a re- 
alistic equation of state, heating and cooling of the disc, and 
a self-consistent treatment of the evolving ionisation fraction 
and conductivity of the disc material. At the present time in- 
clusion of these physical processes is beyond current computa- 
tional resources. 
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